% The code is based on the following paper:
%   Chen, T. Y., Lin, Y. L., and Tzeng, L. Y., forthcoming.
%   Estimating probability weighting functions through option pricing bounds. Review of Asset Pricing Studies.
%
% Copyright: Tzu-Ying Chen, Yo-Lan Lin, Larry Y. Tzeng
% Date: January 30, 2024

clear

%% Setting
RootName_V_COND = ['ATMIV'];
% RootName_V_COND = ['EGARCHV'];

%% Load Data
FileName = ['GrossRET_1981_2021_DiscreteState'];
load(FileName, ...
     'Date_Monthly', 'State_Monthly', 'State_PROB')
clear FileName

%% Volatility Information Used in Conditional Distribution
if strcmp(RootName_V_COND, 'ATMIV')==1
    FileName = ['ATMIV19962021_Monthly_Pre1Day'];
    load(FileName, ...
         'Target_AllDate', 'ATMIV');
    clear FileName 

    [~, Index_Location] = ismember(Date_Monthly(:, 2), Target_AllDate(:, 2));   
    V_COND = ATMIV(Index_Location);                            
    clear Index_Location Target_AllDate ATMIV
    
    AllData = [];                                                          
    for d = 1:size(Date_Monthly, 1)
        Target_Date = Date_Monthly(d, 1);
        Target_TTM = Date_Monthly(d, 3);

        FileName = ['OP' num2str(Target_Date) '_TTM' num2str(Target_TTM) '.mat'];
        load(FileName, ...
             'Data', 'Index_Date', 'Index_TTM');
        clear FileName 

        AllData = [AllData; Data];
        clear Target_Date Target_TTM Data
    end
    clear d

    TTM = unique(AllData(: , [Index_Date Index_TTM]), 'rows');  
    TTM = sortrows(TTM, 1);                                                
    TTM = TTM(:, end);                                                     
    V_COND = V_COND - 0.02;                                                
    V_COND = sqrt(TTM / 365) .* V_COND;                                    
    clear Index_Date Index_TTM AllData TTM
else
end

if strcmp(RootName_V_COND, 'EGARCHV')==1
    V_COND = nan(size(Date_Monthly, 1), 1);                                

    for d = 1:size(Date_Monthly, 1)
        Target_Date = Date_Monthly(d, 1);
        Target_TTM = Date_Monthly(d, 3);    

        FileName = ['DG_GrossRET_Monthly_' num2str(Target_Date)];
        load(FileName, ...
             'RET_Monthly');
        clear FileName

        [AllRET, ~, Index_AllRET] = unique(RET_Monthly); 
        AllRET_PROP = accumarray(Index_AllRET, 1, size(AllRET), @sum, nan) / length(RET_Monthly);
        M1 = sum(AllRET_PROP .* power(AllRET, 1));                         
        M2 = sum(AllRET_PROP .* power(AllRET - M1, 2));                       
        V_COND(d, 1) = sqrt(M2);
        clear Index_AllRET Target_Date Target_TTM RET_Monthly AllRET AllRET_PROP M1 M2
    end
    clear d
else
end

%% Conditional Distribution
for d = 1:size(Date_Monthly, 1)
    State = State_Monthly(d, :);
    M1 = sum(State_PROB .* State);
    M2 = sum(State_PROB .* power(State - M1, 2));                          
    State_Monthly(d, :) = ((State - M1) / sqrt(M2)) * V_COND(d, 1) + M1;   
    clear State M1 M2
end
clear V_COND
clear d

%% Output
Data = State_Monthly;
clear State_Monthly

for d = 1:size(Date_Monthly, 1)    
    State_Monthly = Data(d, :);
    
    Target_Date = Date_Monthly(d, 1);
    FileName = ['CJP_' RootName_V_COND '_GrossRET_Monthly_' num2str(Target_Date)];
    save(FileName, ...
         'State_Monthly', 'State_PROB')
    clear FileName   

    clear Target_Date State_Monthly
end
clear Data State_PROB
clear d
